# Code found and adapted from:
# https://stackoverflow.com/questions/27164114/show-confidence-limits-and-prediction-limits-in-scatter-plot
# Helper functions -------------------------------------------------------------
def equation(a, b):
"""Return a 1D polynomial."""
return np.polyval(a, b)
def plot_ci_manual(t, s_err, n, x, x2, y2, ax=None):
"""Return an axes of confidence bands using a simple approach.
Notes
-----
.. math:: \left| \: \hat{\mu}_{y|x0} - \mu_{y|x0} \: \right| \; \leq \; T_{n-2}^{.975} \; \hat{\sigma} \; \sqrt{\frac{1}{n}+\frac{(x_0-\bar{x})^2}{\sum_{i=1}^n{(x_i-\bar{x})^2}}}
.. math:: \hat{\sigma} = \sqrt{\sum_{i=1}^n{\frac{(y_i-\hat{y})^2}{n-2}}}
References
----------
.. [1] M. Duarte. "Curve fitting," Jupyter Notebook.
http://nbviewer.ipython.org/github/demotu/BMC/blob/master/notebooks/CurveFitting.ipynb
"""
if ax is None:
ax = plt.gca()
ci = t * s_err * np.sqrt(1/n + (x2 - np.mean(x))**2 / np.sum((x - np.mean(x))**2))
ax.fill_between(x2, y2 + ci, y2 - ci, color=tick_color, edgecolor=None)
return ax
# Computations ----------------------------------------------------------------
x = reviews['lengths']
y = reviews['std_rating_hat']
p, cov = np.polyfit(x, y, 1, cov=True) # parameters and covariance from of the fit of 1-D polynom.
y_model = equation(p, x) # model using the fit parameters; NOTE: parameters here are coefficients
# Statistics
n = reviews['std_rating_hat'].size # number of observations
m = p.size # number of parameters
dof = n - m # degrees of freedom
t = scipy.stats.t.ppf(0.95, n - m) # used for CI and PI bands
# Estimates of Error in Data/Model
resid = y - y_model
chi2 = np.sum((resid / y_model)**2) # chi-squared; estimates error in data
chi2_red = chi2 / dof # reduced chi-squared; measures goodness of fit
s_err = np.sqrt(np.sum(resid**2) / dof) # standard deviation of the error
# Plotting --------------------------------------------------------------------
# Plot globals
title_color, tick_color, edge_color, data_color = 'black', 'darkgray', 'dimgray', 'whitesmoke'
font_size = 25
font_weight = "bold"
fig, ax = plt.subplots()
# Data
ax.plot(x, y, "o",
color=data_color,
markersize=.5,
markeredgewidth=.5,
markeredgecolor=edge_color,
markerfacecolor=data_color
)
# Linear Fit
ax.plot(x, y_model, "-",
color=title_color,
linewidth=1,
label="Fit")
# Confidence Interval
x2 = np.linspace(np.min(x), np.max(x), 100)
y2 = equation(p, x2)
plot_ci_manual(t, s_err, n, x, x2, y2, ax=ax)
# # Prediction Interval
pi = t * s_err * np.sqrt(1 + 1/n + (x2 - np.mean(x))**2 / np.sum((x - np.mean(x))**2))
ax.fill_between(x2, y2 + pi, y2 - pi,
color="None",
linestyle="--",
linewidth=3)
ax.plot(x2, y2 - pi, "--",
color=title_color,
label="95% Pred. Lim.")
ax.plot(x2, y2 + pi, "--",
color=title_color)
# Custom legend
handles, labels = ax.get_legend_handles_labels()
display = (0, 1)
anyArtist = plt.Line2D((0, 1), (1, 0), color=tick_color) # create custom artists
legend = plt.legend(
[handle for i, handle in enumerate(handles) if i in display] + [anyArtist],
[label for i, label in enumerate(labels) if i in display] + ["95% CI"],
loc=9, bbox_to_anchor=(1, .6, .5, 0.), ncol=1, mode="expand", fontsize=20)
frame = legend.get_frame().set_edgecolor(tick_color)
# Annotation with r
plt.text(2200, 2.4, 'r = {}'.format(round(corr[0],2)), fontsize=25)
# Styling
plt.xlabel('N Tokens per Review',
fontsize=font_size,
labelpad=25,
color=title_color,
weight=font_weight)
plt.ylabel('Predicted\nRating',
fontsize=font_size,
rotation=0,
labelpad=90,
color=title_color,
weight=font_weight)
plt.grid(b=True,
color=tick_color,
axis="y",
alpha=0.8,
linestyle=':',
linewidth=1)
plt.xticks(fontsize=25,
color=tick_color,
weight=font_weight,
rotation = 60)
plt.yticks(fontsize=25,
color=tick_color,
weight=font_weight)
plt.xlim(0, 3000)
plt.ylim(-3, 3)
plt.show()
